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We used a quantum Monte Carlo method to study the magnetic impurity adatoms on graphene. 
We found that by tuning the chemical potential we could switch the values of the impurity's local 
magnet moment between relatively large and small values. Our computations of the impurity's 
spectral density found its behavior to differ significantly from that of an impurity in a normal metal 
and our computations of the charge-charge and spin-spin correlations between the impurity and 
the conduction band electrons found them to be strongly suppressed. In general our results are 
consistent with those from poor man's scaling and numerical renormalization group methods. 
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I. INTRODUCTION 

Graphene is a two-dimensional fermionic material 
whose band structure has a pseudo-gap created by a 
particular arrangement of touching Dirac cones.— In the 
vicinity of these cones, the electronic density of states 
p(E) varies linearly with the energy E measured relative 
to the Fermi energy Ep, that is, p(E) = a.\\E — Ep\. 
This functional behavior and the corresponding low den- 
sity of states open graphene to the possibility of tailoring 
unconventional behavior. It is well known, for example, 
that magnetic impurities behave differently in a pseudo- 
gapped material than in a normal metal.— ^ 

In a normal metal, the magnetic impurities induce 
many-body correlations that at low temperatures quench 
the spin fluctuations at the impurity site. This is the 
Kondo effect. For this phenomena two impurity mod- 
els, the Anderson and Kondo models, have been partic- 
ularly well studied by poor man's scaling and numerical 
renormalization group methods^ Both methods identify 
fixed point Hamiltonians of the same form as the origi- 
nal Hamiltonians but with renormalized parameters. For 
the Anderson model the original parameters are £<j, U, 
and r, which are the impurity's energy level, Coulomb 
interaction between two electrons simultaneously occu- 
pying this level, and level width. The renormalization 
mainly affects £a, and all renormalizations flow to the 
only stable fixed point, the strong coupling fixed point, 
accompanied by Tximp — > and Ximp — > constant as the 
temperature T — > 0. The thermodynamic and transport 
properties of flows passing in the vicinity of the unstable 
local moment fixed point exhibit universality when their 
temperature and frequency dependencies are scaled by 
the Kondo temperature Tk- 

The eigenvalues and eigenvectors of strong coupling 
limit of the Anderson model, when only one electron oc- 
cupies the impurity level, are equivalent to those of the 
Kondo model in its weak exchange limit.— The equiv- 
alence establishes a well-defined relation between the 



effective Anderson and the actual Kondo exchange in- 
teractions between the impurity moment and the con- 
duction electrons. For the Kondo model with an anti- 
fcrromgnetic exchange J > 0, the renormalizations al- 
ways flow from its unstable J = fixed point to its stable 
J = oo fixed point at which its local moment is quenched. 

Work on both models has typically assumed that 
the conduction band density of states p(E) is constant. 
There however is now a well established body of litera- 
ture that shows when p(E) = ev r |-E| r (setting Ep = 0) 
and r > 0, the two models have features with no counter- 
parts in the constant p(E) models.— ££~— For the Kondo 
model, the J = fixed point becomes stable and a new 
unstable J c fixed point appears. The Kondo effect oc- 
curs only if J > J c but J c can be very large when r > , ■ 
For the Anderson model, as r increases, the effective ex- 
change interaction weakens and very difficult to boast 
above J c . The once stable strong coupling fixed point 
becomes unstable and the the local moment fixed point 
becomes stable so at T = partially quenched moments 
can survive. Further both Ed and Y are renormalized, and 
both spin and charge fluctuations are suppressed. Also 
the stability and nature of some fixed points depend on 
whether particle-hole symmetry exists. Thus, for pscudo- 
gapped materials Kondo quenching of the magnet impu- 
rity often will not exist. 

Clearly, the Kondo difficult case where r = 1 , for which 
most of the just stated phenomena occurs, is relevant 
to graphene. Indeed, several studies exist that focus 
on graphene as an opportunity to study the Kondo ef- 
fect and Kondo quantum criticality in a pseudo-gapped 
material Air— There is also considerable interest in ex- 
ploiting the now well-established experimental capability 
of shifting graphene's chemical potential p by an applied 
electric field to switch on and off this novel physics J2r— 

In this paper we share the interest in using an electric 
field to switch the properties of graphene. To this end 
we studied the Anderson impurity model for graphene as 
a function of T and p by using a determinant quantum 



Monte Carlo method based on Hirsch-Fye algorithm^ 
Instead of using p(E) = ai\E\, we used the actual density 
of states for a tight-binding expression of graphene's con- 
duction band. We note that when p ^ 0, p(E) is replaced 
by p(E — p), destroying the symmetry p(E) = p{~E) as- 
sumed by scaling and renormatization methods. Hence 
besides exploring the cases with variations from linearity, 
we are also exploring cases where the conduction band 
density of states is not symmetric about Ep . 

Wc find that over a reasonably wide range of parame- 
ters a local moment, in the sense of a non-zero expecta- 
tion value of S%, exists. As Ep is gated to below zero, 
the renormalized impurity level e* d simultaneously shifts 
toward it. Eventually, the two energies pass each other, 
transferring charge from the impurity to the conduction 
band and in the process decreasing the magnitude of the 
moment on the impurity. The process thus "switches" 
the magnetic moment from a relatively large value to a 
relatively small one as a function of the gating. In fact 
the switch is from a relatively well-developed moment 
to one that is partially screened. Our computations of 
the spectral density of the impurity support not only the 
shifting of e* d but also a significant reduction in the value 
of r. These changes are consistent with the renormal- 
ization group's results of reference.— We note however 
the calculations there were only for the p = case. We 
expect the p needs to be moved out of the linear density 
of states region before we would see the renormalization 
of r to cease. We also compute the charge-charge and 
spin-spin correlations between the impurity and conduc- 
tion band electrons and find them to be small amplitudcd 
and short ranged. 



II. FORMULATION 

The Anderson impurity model with single impurity or- 
bital of energy Ed and Coulomb repulsion U, couples the 
conduction electron states and impurity with hybridiza- 
tion V. The total Hamiltonian is H = Hq + Hi + E 2 - 
Hq is a tight-binding Hamiltonian. For graphene it is 



Here (ij. creates an electron with spin a at the impurity 
orbital. Finally H 2 describes the hybridization between 
the impurity adatom and a graphene atom 

H 2 = Vj2iald a +4a 0a }. 

a 

We assume the impurity is on the top of the site Ro a of 
sub-lattice A. 

As previously noted, wc simulated this model 
with the Hirsch-Fye quantum Monte Carlo algorithm. 
This algorithm^! naturally returns the imaginary-time 
Green's function Gci(t) = ^Gda (t) of the impurity. 

(7 

With this Green's function we determined its associated 
spectral density A(u) = ^ CT A a (uj) by numerically solv- 
ing 



^S^T- 



Specifically, we used the procedures detailed in Ref.l23lfor 
qualifying the data and qualifying the solution. We found 
that the three different Bayesian methods for doing the 
analytic continuation described in Ref. [23| produced only 
small differences in computed A(uj) and that the results 
were similarly insensitive both to the use of Gaussian 
and flat default models and to the use of a constant and 
Jeffery prior. 

We also used an extension of the Hirsch-Fye 
algorithm^ to compute the charge-charge correlation 
function 

Ci = (n d n. l ) - (n d )(ni), 

and the spin-spin correlation function 

Si = (m d mi}, 

where rii and rrii are the charge and magnetic moment of 
the graphene atom at site i. 



H = -t J2 [alb jc 



AW] 



vYJ a \<? a i<y + b Lh 



<n>,° 



where a itT and b ia creates an electron with spin a at sites 
Ri a and R^ on the A and B sub-lattices of graphene's 
hexagonal structure. In graphene the hopping matrix 
clement t > is about 2.8 eV— and p is chemical potential 
to be tuned by a gate voltage. There arc two bands, the 
■k and 7r* bands, each of width of 3£, that touch each 
other at six Dirac points in the first Brillouin zone of a 
hexagonal lattice. When p = 0, the density of states near 
E F = is p{E) = ai\E\ with ai = 4V3/37rf 2 . H x is the 
impurity Hamiltonian 



H X = 5> d - p)4d a 



Udld^dld^ . 



III. RESULTS 

A. Magnetic Moments 

In Fig. QJi-c, we show various physical quantities as a 
function of p for different values of V. These are the 
impurity level occupancy rid = (n^ + ridi), its double 
occupancy n np n ( j L own = {nd^nd\), and its local moment 
squared m d = ((ndf — ndx) 2 ). To the accuracy of our sim- 
ulation m d = (ridt - n di) = 0, implying {n d -\) = (ridi)- 
All three quantities "switch" their values as the chemical 
potential moves below the Dirac point. For p near this 
point, each case tends to an average occupancy of 1 but 
even for this case we note that some double occupancy 
is present. Also noting that m d = n d — 2n up ndown, we 
see that the switching of m d is driven by the switching 
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FIG. 1: (Color online), (a) Occupancy rid, (b) double oc- 
cupancy nd^ridi, (c) local moment squared m d , and (d) the 
susceptibility T\ versus the chemical potential /i. V is the 
hybridization, Ed = — U/2 — — 0.40£, and the inverse temper- 
ature T" 1 = 64t -1 . 



in rirf accompanied by some reduction in n up «down. It is 
interesting that for different values of V the tunable re- 
gions occur over different ranges of ^ and that the smaller 
values of V produce the bigger effects but require larger 
values of /x. Below we will connect much of this behavior 
with e* d and T* shifts as V changes. 

We also calculated the temperature dependent, impu- 
rity spin susceptibility^ 



X(T) 



f 

Jo 



dr(m d (r)m d (0)), 



where /3 = T" 1 and m d (T) = e TH m d (0) e - TH . In Fig.EEJl 



we show Tx versus \x. Clearly, its behavior correlates 
with that of m d . Figure [5] shows \ as a function of T _1 
for various values of /z. Here V = l.Oi, T _1 = 64i _1 , 
U = 0.80i, and Ed = — 0.40i. As /j, moves below the Dirac 
point and T is lowered, we see that x crosses over from a 
Curie- Weiss behavior to the behavior of an screened local 
moment. 
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FIG. 2: (Color online) The spin susceptibility \ versus the in- 
verse temperature T _1 for various values of chemical potential 
fi. In all cases V = l.Oi, U = 0.80i, and e d = -0.40t 



B. Spectral Densities 

In Fig. [3] are the spectral densities A(uj) for an inverse 
temperature 1/T = 12t _1 and e d = -U/2. In Fig. [3Ji 
we fix /i at and V^ at l.Oi and vary U. At /j, = 0, the 
symmetry of the bands and the choice £d = —U/2 places 
the Anderson model in a state of particle-hole symmetry. 
This symmetry implies A{uj) = A(— uj) which is evident. 
The pseudo-gap is also evident. Additionally, we see that 
as U increases the two peaks of A(ui) increase their sepa- 
ration and broaden. As the peaks broaden, their heights 
collapse to accommodate the sum rule / A(uj) du; = rid. 

The- features of the A(oj) in Fig. [3] differ markedly from 
several general features of a Hartree-Fock solution for a 
normal metal^ where the peak heights and widths are 
controlled by V and independent of U and their separa- 
tion D ss U. The exact results on the other hand has 
peak heights and widths varying with U and peak sepa- 
rations I? at a given value of U being much smaller than 
U. 

Figure [3}d shows A(u) for different hybridizations V 
but with /i still equal to and Ed still equal to —U/2. 
Here we see additional differences from Hartree-Fock for 
a normal metal: When V increases, the A(uj) peaks shift 
toward the Dirac point and become sharper and higher. 
This behavior is consistent with Hartree-Fock calcula- 
tions using a linear density of states but oppositely trends 
the predicted behavior of Hartree-Fock calculations with 
a constant density of states where increasing V makes the 
peaks broader and lower. Additionally, for the symmetric 
model, the peak positions do not shift. 



6 


(a) 


» 


1* 


j 


L ■ 


2 


i 


i 


n- 


M 


UIv . 




-0.4 -0.2 0.0 0.2 0.4 
to/t 



|6 
< „ 





-0.4 


-0.2 


0.0 
co/t 


3.2 


0.4 






• 


A 


A 










0.6' 


■ 
• 
▲ 


U 
U 

u 


=1.6t 
=1.2t 

=0.80t 


. 


3 

0.4. 


- 




A 


. 


0.2. 


■ 


• 
■ 


• 
■ 


• 
■ 




A 

• 

• 

■ ■ 


on. 















-0.30-0.15 0.00 0.15 0.30 0.45 0.60 0.4 
to/t 



FIG. 3: (Color online), (a) The spectral density A(u) versus 
us for V = l.Ot, /i = 0, and (from top to bottom) U = 0.80i, 
1.2i, 1.6*, and 2.0i. (b) A(ui) versus w for U = 0.80i, /i = 0, 
and (from top to bottom) V = l.Oi, 0.75i, 0.60i,and 0.50i. 
(c) A{uj) versus u> for jj, = — 0.15i, U = 0.80i, and (from top 
to bottom) V = 0.75i, 0.60i, and 0.50i. (d) The distance 
D between two peaks of A(u>) versus V/U and fi — 0. Here 
T' 1 = 12i _1 . In all cases e d = -f//2. 
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FIG. 4: (Color online). The charge-charge d and spin-spin 
Si correlations versus site i. V = 0.75*, T -1 = 64£, and 
£d = —U/2 — — 0.40£. The adatom is at i — and the 
numbering of the other sites is given in (a). The insets of d 
and Si show the details for tails of curves. One inset to (b) 
shows the spin-spin correlation in real space with fi = 0. The 
black ball represents the impurity adatom at an A sub-lattice, 
and the filled red and open blue circles represent lattice sites 
with negative and positive values of Si at a distance i from 
the black ball. 



In Fig. [3]: we examine the case of gated impurity-doped 
graphenc; that is, we eliminate the particle-hole symme- 
try in p(E) by having li = — 0.15*. We see that the asym- 
metric A(ui) is enhanced when w is negative, and both 
peaks display trends similar to those found in Fig.[5}D but 
the uj < peaks are sharper, higher, and closer to u> = 0. 

Finally, in Fig.[3Jl, we summarize the energy difference 
D in detail. For a fixed Ed = —U/2, increasing V de- 
creases D. Fixing V and increasing U{= —2ed) decreases 
D. For example, when U = 0.80t and V = 0.40C7, D is 
about 70% of U while U = 1.6* and D is only about 28%. 
The bare Ed differs so much from its renormalized value 
e* d that the impurity level may be detectable even if the 
Ed seems outside the experimentally accessible range. 

We can loosely correlate the li and V dependences of 
the A(cj) in Fig. [3J: with those of rid, rn d , and n up Kdown 
in Fig. [TJ In Fig. [T] we see that switching occurs well 
before Li reaches Ed = — 0.40*. We also see that the 

V = 0.50*, 0.60*, and 0.75* cases at li = -0.15* corre- 
sponds to unswitched, just started switching, and switch- 
ing cases. Comparing Fig. 0J: with Fig. Q] reveals that at 

V = 0.75* starts switching just as li is dropped past the 
A(u) peak at w « -0.10*. At V = 0.60*, /x = -0.15* is 
approximately the value of the frequency at the left edge 
of A{uj), and &tV = 0.50*, li = —0.15* sits at the peak. 



C. Spin and Charge Correlations 

The linear energy dispersion and the vanishing of the 
density of states at the Dirac points generate for impuri- 
ties an unusual Friedel sum rulej^i Friedel oscillations^ 
and RKKY interaction<22r— The correlations of the im- 
purity spin and charge with those of the conduction elec- 
trons reflect these behaviors. For example, when both 
H = and U = 0, instead of a Fermi surface, graphene 
has two Fermi points at the two non-equivalent Dirac 
points (that is, the Dirac points not connected by a re- 
ciprocal lattice vector). Perfect nesting exists between 
these pointSp2r— leading to spin and charge densities 
without oscillations. The magnitude of the nesting wave- 
vector is K = 47r/3v / 3a, where a is the carbon-carbon 
spacing. Predicted for RKKY interactions, for example, 
are short-ranged ferromagnetic correlations between the 
local moment and the conduction electron spins instead 
of the standard anti-ferromagnetic one and an oscillation 
pattern, determined by K, such that if the impurity is at 
an A sub-lattice site and so is i, the sign of these oscil- 
lations is negative, and if i is at a B sub-lattice site, the 
sign is positive. 

In Fig. @J for V = 0.75*, we present examples of the 
behavior of Cj and Si when U ^ for cases when /i is 
zero and not zero. In these figures the impurity adatom is 
located on the top of the site labelled 0. The subsequent 
labeling of the lattice sites is shown in the inset to Fig. 4a. 
When [x = 0, we see that the charge correlations still lack 
oscillations, but the formation of a local magnetic mo- 
ment (Fig. [T]) leads to oscillating spin corrclation o 29 ' 31 ' 33 
on a length scale set by K. The nearest neighbor spin 
correlations are ferromagnetic instead of the standard 
anti-ferromagnetic correlation.— When li ^ 0, both the 
spin and charge correlation functions still appear to os- 
cillate on a scale set by K instead of twice the Fermi 
wave- number kp where kp is defined by \li\ = vpkp 
where vp = 3*/2a is the Fermi velocity. Predications 
have this Fermi scale interfering^ or dominating 3 -^ the K 
scale. Fig. [4] shows that the short-ranged correlations re- 
vert to the standard anti-ferromagnetic ones. The ll 7^ 
oscillation pattern appears phase shifted relative to the 
particle-hole symmetric case. We lack accuracy to iden- 
tify interference with a 2kp scale. 

In general, for the interacting problem the spatial ex- 
tent of both the spin and charge correlations is relatively 
short ranged, and their amplitudes are small. Doping 
most clearly changes the correlations in close proxim- 
ity to the impurity. When li is in the region of linear 
electronic dispersion, the length scale of the oscillations 
reflects the geometric length scale and not the doping. 

We can also loosely correlate the li dependence of the 
these correlation functions with that of rid, m d , and 
iup«down in Fig. [1] If we follow the ll dependence of 
the V = 0.75* curve in Fig. QJt and lc, we see that 
near \x = all three quantities have their maximum val- 
ues. When li is near —0.20* and —0.30*, all three val- 
ues drop. Overall the drops in rid an d m d lead to de- 



creased correlations. The decrease in n U p«down creates 
a stronger on-site Fermi-hole effect and hence stronger 
on-site anti-correlations. The drops between fj, = — 0.20t 
and fj, = — 0.30i are relatively small and hence only create 
small changes in the correlation functions at these values 
of /i. 



IV. CONCLUSIONS 



is only partially screened. We found unconventional be- 
havior for the impurity spectral densities and correlation 
functions that further highlight the difference between an 
impurity in pseudo-gapped graphene and one in a metal. 
We suggest that a scanning tunneling microscopy (STM) 
can measure the spectral densities and the charge-charge 
correlation functions and a spin-polarized STM can mea- 
sure the spin-spin correlations^ - — 



In summary, our calculations support prior suggestions 
that it should be possible to switch the magnetic moment 
of an impurity adatom on the surface of graphene from 
a relatively high value to a relatively low one by shifting 
the chemical potential by an electric field. Being shifted 
is a reasonably well-defined local moment to one that 
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